Analysis of directed flow from elliptic flow 
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The directed flow of particles produced in ultrarelativistic heavy ion collisions at SPS and RHIC 
is so small that currently available methods of analysis are at the border of applicability. Standard 
two-particle and flow-vector methods are biased by large nonflow correlations. On the other hand, 
cumulants of four-particle correlations, which are free from this bias, are plagued by large statistical 
errors. Here, we present a new method based on three-particle correlations, which uses the property 
that elliptic flow is large at these energies. This method may also be useful at intermediate energies, 
near the balance energy where directed flow vanishes. 

PACS numbers: 25.75.Ld 25.75.Gz 



I. INTRODUCTION 

The azimuthal angles of particles emitted in a heavy 
ion collisions are correlated with the direction of the im- 
pact parameter (reaction plane) of the two incoming nu- 
clei. This correlation is characterized by the Fourier co- 
efficients of the particle azimuthal distribution, v n B: 



v n = (cos(n(0 - 



(1) 



where n is an integer, <f> denotes the azimuthal angle 
of an outgoing particle, (5>r is the azimuth of the reac- 
tion plane, and angular brackets denote an average over 
events and over all particles in a given transverse mo- 
mentum and rapidity window. 

The first Fourier coefficient V\, referred to as directed 
flow, is usually positive in the forward hemisphere. This 
is because particles are deflected away from the target. 
This effect was first observed at Bevalac in 1984 0] , and 
reaches its maximum value at about 400 MeV per nucleon 
||. At much lower energies, below 100 MeV per nucleon, 
v\ becomes negative due to the attractive nuclear mean 
field j|, |f|. At ultrarelativistic energies, it decreases: at 
the CERN Super Proton Synchrotron (SPS), it is of the 
order of 2% g. At the Brookhaven Relativistic Heavy 
Ion Collider (RHIC), it has not been measured yet. Ac- 
curate measurements of V\ are important since this quan- 
tity is a sensitive probe of nuclear matter properties |7j. 
There is an interesting prediction that this quantity may 
even become negative above SPS energies ||. 

Standard methods for analyzing directed flow |], [l(], 
[HJ are based on the assumption that all azimuthal corre- 
lations between outgoing particles are due to flow. When 
flow becomes too small, however, other sources of correla- 
tion ("nonflow" correlations) can no longer be neglected. 
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Correlations due to global momentum conservation are 
well known jL2j, and a systematic method to correct for 
this effect has recently been proposed However, 
many other effects are expected to produce nonflow corre- 
lations of the same order of magnitude, such as quantum 
(HBT) correlations jljj and resonance decays Jl^] . We 
have proposed new methods, based on higher order cor- 
relations jl6[ |lj] , which allow one to get rid of nonflow 
correlations. However, statistical errors on four-particle 
(or higher order) correlations are larger, so that these 
methods may fail if flow is too small. 

In this paper, we propose a new method which is free 
from these limitations. It allows one to obtain an esti- 
mate of i>i which is insensitive to nonflow correlations, 
with statistical errors usually not much larger than with 
standard techniques. This method is based on the ob- 
servation that elliptic flow (v2, second Fourier harmonic 
of the azimuthal distribution) is reasonably large in the 
situations considered where directed flow is hard to mea- 
sure: V2 is roughly 3% at SPS p], and even larger at 
RHIC [Jis| . Our method relies on the measurement of 
three-particle azimuthal correlations, which involve both 
vi and V2- 



(vi) 2 v 2 , 



(2) 
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where 4>i, <f>2, and ^3 denote azimuthal angles of three 
particles belonging to the same event, and the average 
runs over triplets of particles emitted in the collision, and 
over events. This equation yields (vi) 2 i>2, and one then 
deduces v± once V2 has been obtained from a separate 
analysis. 

The idea of mixing two Fourier harmonics is not new. 
It underlies measurements of elliptic flow relative to the 
event plane determined from directed flow p9[ : in these 
analyses, one first obtains v\, and then V2 from an equa- 
tion similar to Eq. (||). Here we propose to do the other 
way round: first measure V2, then obtain v\ from the 
three-particle correlation. 

In Sec. H 

we present in more detail the principle of 
the method, and we explain why it is better than other 
methods when directed flow is small. Then, in Sec. Ill, 
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we describe the practical implementation of the method. 
This implementation makes use of two techniques already 
introduced in Refs. [ p"7| p0| |: first, a cumulant expansion 
allows us to eliminate the effects of azimuthal asymme- 
tries in the detector acceptance; second, the formalism 
of generating functions provides an elegant way of con- 
structing three-particle correlations (and beyond), with a 
computer time which grows only linearly with the num- 
ber of particles involved, while taking into account all 
possible particle triplets. A generalization to higher or- 
der cumulants is presented in Sec. IV. Our results are 
discussed in Sec. [v|. 



II. PRINCIPLE AND ORDERS OF 
MAGNITUDE 

In this section, we compare three methods of analyz- 
ing directed flow: the standard two-particle technique, 
the four-particle cumulant, and the three-particle mixed 
corre latio n, Eq. (^), which is the focus of this paper. In 
Sec. II A , we recall how vi can be obtained from 

Sec " 



It is straightforward to generalize Eq. (|3|) to higher 
order correlations, such as the four-particle correlation: 
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imuthal correlations in various ways. In 
Sec. II C, we estimate for each method the order of mag- 
nitude of errors due to nonflow correlations and finite 
statistics. Numerical estimates are given in Sec. [ID , 
where we show that three-particle correlations provide 
the most (if not only) reliable way of an alyzin g directed 
flow at ultrarelativistic energies. In Sec. HE , we finally 
explain how to perform detailed measurement of V\ as a 
function of transverse momentum and rapidity with this 
method. 



A. Directed flow from azimuthal correlations 

In a given collision, the azimuth of the reaction plane 
is unknown, so that the Fourier coefficients v n of the 
azimuthal distribution, defined in Eq. (|l|), can only be 
obtained from azimuthal correlations between the out- 
going particles. The simplest way is to use correlations 
between two particles |pX| , labelled 1 and 2, belonging to 
the same event : 



* , l-*R)gi(*R-02) 



3 i(01-$ H ) 



3 i(* R -0 2 ) 



= M 2 , (3) 

where, as in the previous equations, averages are taken 
over pairs of particles and over events. In going from 
the first to the second line, we have assumed that all az- 
imuthal correlations are due to flow, i.e., that azimuthal 
angles relative to the reaction plane <j)\ — $_r and 4>2 — 
are independent. The "standard" flow analysis [|| pro- 
ceeds differently: one correlates the azimuth of a parti- 
cle with a "flow vector" obtained by summing over many 
particles: this involves a sum of two-particle correlations, 
and this method is in fact essentially equivalent to the 
two-particle technique. 



A't>l+<t>2-<t>3-<t>i)\ _ ( Wl ) 4 ; 



(4) 



where the average now involves all possible 4-uplets of 
particles belonging to the same event. This result will be 
used below. 

In this paper, we shall be concerned with yet another 
type of correlation, namely, a three-particle correlation 
which mixes the first two Fourier harmonics: 



e i(<t>l+4>2-2<p3)\ _ / e i(0l-*H) e i(<fe-*H) e 2i(*B-03) 



(vi) v 2 . 



(5) 



If V2 is already known from a previous analysis, using the 
three-particle correlation Eq. (||) allows one to extract an 
estimate of v\. 



B. Nonflow correlations 

The above estimates were derived under the assump- 
tion that all azimuthal correlations are due to flow. How- 
ever, there arc also other, "nonflow," contributions to 
azimuthal correlations due to various effects 0, H, [|. 
Let us recall the order of magnitude of these correlations 
by taking a simple example: if a resonance decays into k 
particles, these k particles will be strongly correlated by 
the decay kinematics. Now, the probability that k arbi- 
trary particles seen in a detector originate from the same 
resonance scales like the total multiplicity M of the event 
like 1/M k ~ 1 . This is generally the order of magnitude of 
the genuine fc-particle correlation due to nonflow effects. 
In particular, the two-particle nonflow correlation is of 
the order of 1/M, so that Eq. (B|) reads in fact 



(6) 



When vi becomes smaller than 1 / V M, the resulting error 
on v\ due to the unknown nonflow term becomes as large 
as vi itself. This is probably the case at SPS, as shown 
in Ref. |lq |, and a fortiori at RHIC energies. 

A similar reasoning applies to higher order correla- 
tions. In the case of the four-particle correlation, Eq. (Q), 
it may for instance happen that two pions labelled 1 and 
3 originate from the same p meson, while 2 and 4 come 
from another p meson. This gives a four-particle correla- 
tion of order 0(1/M) 2 . However, these pairwise correla- 
tions can be subtracted from the measured four-particle 
correlation so as to isolate the genuine four-particle cor- 
relation. This is the principle of the cumulant expansion 
which was proposed in Refs. jl6, 17 to get rid of nonflow 
correlations in the flow analysis. The cumulant of the 
four-particle correlation is defined as 
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i(01+02-03-04) \\ — / „i(01+02-03-</>4) 



= i(01-0 3 ) 



5 i(02-04) 



= i(0i-0 4 ) 



,i(0 2 -03; 



(7) 



r 



According to the above discussion, the contribution of 
nonflow effects to this genuine four-particle correlation 
is of order l/M 3 , much smaller than 1/M 2 . The con- 
tribution of flow follows in a straightforward way from 
Eqs. (H) and (0), so that one may finally write 



| e i(*i+*>-*>-#*)| = _( Wl )4 + 



(8) 



Finally, let us estimate nonflow contributions to the 
mixed three-particle correlation, Eq. (|5|). Unlike the 
four-particle correlation Eq. (^J), this quantity does not 
receive any contribution from two-particle correlations, 



since quantities such as (e l 
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) 



ish by symmetry. The only nonflow correlation is the 
genuine three-particle correlation, of order 1/M 2 , and 
Eq. (||) becomes 



D i(0l+02-2</. 3 )\ _ 



(9) 



In the following, we shall denote by «i{2}, wi{4| and 
i>i{3} the estimates of vi obtained from Eqs. @, m) and 
(^) , respectively, ignoring the nonflow term. Using these 
equations, one finds that the differences due to nonflow 
correlations between these estimates and the exact value 
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vi are of order: 






Vl {2} 


- Vi 


= O 


«i{4} 


- Vl 


= O 


M3} 


- Vl 


= O 



Mvi 
1 



(M Vl ) 3 
1 



{Mv l ){Mv 2 ) 



(10) 



As explained in Ref. (T(| , it is possible to measure v n only 
if Mv n 3> 1, and we assume throughout this paper that 
this condition holds both for v\ and i>2- Then, Eq. ( |l0| ) 
shows that estimates of v\ from three- or four-particle 
correlations are much less biased by nonflow correlations 
than standard estimates from two-particle correlations. 



C. Statistical errors 

In practice, the use of higher order correlations is lim- 
ited by statistical errors due to limited statistics. The 
various correlations encountered are quantities of the 
type (1/iV) Ylj=i cos(A^j), where the A(j>j are various 
combinations of the particle azimuths. If the A(j)j are in- 
dependent and randomly distributed, the standard error 
on such a quantity is l/\/2N. 



In the case of two-particle correlations, Eq. (|3|), one 
can construct M(M — l)/2 different pairs of particles in 
an event with multiplicity M, so that the total number 
of combinations is N ~ N cvts M 2 /2 for 7V evt s events. In 
the case of four-particle correlations, Eq. (||), one can 
construct M(M - 1){M - 2)(M - 3)/8 independent 4- 
uplets, so that N ~ iV cvts M 4 /8. Finally, in the case of 
three-particle correlations, there are M(M — 1)(M — 2)/2 
independent triplets, so that N ~ N cvts M 3 /2. 

One thus obtains the following expressions for the rel- 
ative statistical errors on the various estimates of v\: 



fai{4} 
fai{3} 

Vl 



1 



1 



1 1 



2VAU; (wiVM) 4 

1 1 



2VAW^ ( Ul% /M) 2 ( V2 VM)' 



(11) 



These estimates are correct as long as Vi\M and vi\ M 
are not larger than unity, otherwise correlations due to 
flow must be taken into account in the calculation of 



statistical errors 17 . A mo re acc urate formula for <5vi{3} 
will be given below in Sec. [II D . 

In this paper, we are interested in the situation where 
directed flow is very small, vi\f~M <C 1, which is the case 
when the standard fl ow an alysis fails due to large nonflow 
correlations (see Sec. [IB ). In that case, one sees that the 
statistical uncertainty on i^^}, i.e., the estimate from 
the fourth-order cumulant, is much larger than the error 
on the standard estimate wi{2}. On the other hand, if el- 
liptic flow is significantly large, i^a/M is not much lower 
than unity: in such a case, the mixed-correlation tech- 
nique provides an estimate of v\ with a statistical error 
of the same order as the standard analysis, but which is 
much less biased by nonflow correlations. 



D. Orders of magnitude at SPS and RHIC 

Let us estimate numerically the various errors dis- 
cussed above in a realistic situation. In practice, eval- 
uating the magnitude of nonflow correlations requires a 
detailed modeling. Equations ( |l0| ) , strictly speaking, rep- 
resent scaling laws rather than orders of magnitude. De- 
tailed studies of various effects at SPS energies |l3|, |lj, 
show that the contribution of nonflow correlations to vi 
obtained from the standard analysis (which is equivalent 
to «i{2} above) is of the same order as v\ itself. For 
mid-central collisions, i>2 is about 3% || and the number 
of detected particles M ~ 300. According to Eq. Jic|), 
one expects that using three-particle correlations (i.e., 
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«i{3}) will reduce systematic errors due to nonflow cor- 
relations by a factor of at least Ms 2 , that is, a factor of 
10. Thus one expects estimates using higher order cor- 
relations «i{3} and wi{4} to be little biased by nonflow 
correlations, contrarily to «i{2}. 

Statistical errors can be estimated quantitatively using 
Eqs. (pd|). Realistic values at SPS are iV cvts = 50k events, 
M = 300 particles, V\ ~ 2% and v 2 — 3%. One then 
obtains 



5vi{2} 

v\ 
foi{4} 

v\ 
fai{3} 

Vi 



2% 

16% 

4%. 



(12) 



The relative statistical error on wi{4} is too large to al- 
low detailed measurements of v\ as a function of pr or 
y. On the other hand, the uncertainty on fi{3} is only 
twice larger than the error on Wi{2}, while we have seen 
above that the gain on the systematic error due to non- 
flow correlations is typically a factor of 10. 

At RHIC, v 2 is larger, typically 5% for mid-central 
collisions [Q, so that measuring directed flow from el- 
liptic flow (wi{3}) is even more appropriate than at SPS. 
This is reflected in both systematic errors due to nonflow 
correlations, and in statistical errors. While v 2 is larger 
than at SPS, v\ is expected to be smaller if one extrap- 
olates the decrease observed at SPS compared to AGS. 
The number of detected particles M can be estimated us- 
ing the values of v 2 and the event plane resolution given 
in Ref. [|l8), and is similar to that used above for SPS, 
M ~ 300. 

Since v% is smaller than at SPS, one expects from 
Eq. @ that the bias on fi{2} from nonflow correla- 
tions will be even worse. However, the decrease in v\ 
may be partially compensated by the increase in v 2 , so 
that the error on fi{3} remains of the same order. With 
Aevts = 50k events, V\ = 1%, v 2 = 5%, one obtains the 
following statistical errors 



^i{4} 

Vl 

Svi{3} 

Vl 



7% 

250% 

11%, 



(13) 



where we have used the formula derived in Sec. HID for 
fai{3} [here, v 2 \f~M is of order unity, so that the third of 
Eqs. ( jll| ) no longer applies.] The statistical uncertainty is 
only 50% larger on wi{3} than on vi{2}, while once again 
the gain on the systematic error more than compensates 
for that loss. 

These numerical estimates show clearly that three- 
particle mixed correlations offer the best compromise to 
measure Vi at ultrarelativistic energies. The correspond- 
ing estimate fi{3} is much less affected by nonflow cor- 



relations than standard two-particle methods. In this re- 
spect, it shares the advantages of higher order estimates 
such as «i{4}. In addition, i>i{3} is much less limited by 
statistics than the latter, which would require millions of 
events at SPS and RHIC energies. 



E. Integrated flow, differential flow 

As with other methods, the analysis proceeds in two 
steps. One first estimates the average value of v\ over 
phase space. This is done by averaging e t (<l>i+4>2- 2 <t>3) over 
all possible triplets, as in Eq. (||). However, v\ and v 2 de- 
pend strongly on rapidity and transverse momentum (for 
instance, v± has opposite signs in the backward and for- 
ward hemispheres). In practice, one therefore performs 
a weighted average, and Eq. (0) becomes 



w 1 (l)w 1 (2)w 2 (3)e l( ^ + ^- 2 ' t ' 3 A = (w lVl ) 2 (w 2 v 2 ) , 



(14) 

where w% and w 2 are weights appropriate to directed flow 
and elliptic flow, respectively, which may be any func- 
tion of the particle type, its transverse momentum pt 
and rapidity y. In this equation, w n (k) is a shorthand 
for w n (pTkiVk)- The right-hand side (rhs) of Eq. jl4| ) 
naturally involves weighted averages (w n v n ), rather than 
v„. 

The best choice for the weights is that which leads to 
the sm allest statistical errors. Repeating the discussion 
in Sec. II C, one e asily shows that this is done by maxi- 
mizing (w n v n ) I \/ {Wn). Therefore the best weight is the 
flow itself |l6|, [H), w n (p T ,y) = v n (p T ,y), where v n (p t ,y) 
denotes the value of the flow in a small (j>t, y) bin. In 
practice, one can choose as a first guess the center-of- 
mass rapidity for directed flow, w\ — y — ycM, and the 
transverse momentum for elliptic flow w 2 = pt, in re- 
gions of phase space covered by the detector acceptance. 

Using the value of the integrated (and weighted) el- 
liptic flow (w 2 v 2 ) obtained from a separate analysis, one 
finally obtains the integrated directed flow (wiVi) from 
Eq. (|l4|). Naturally, one must use the same weight w 2 
in the reference analysis which gives (w 2 v 2 ) and in the 
mixed-correlation analysis. 

The second step is to analyze differential flow, i.e., to 
obtain values of v± as a function of transverse momentum 
Pt and/or rapidity y. For that purpose, one averages 
e i(0i+02-20 3 ) over a u 2 an( j b ut res tricts fa to a 

given particle type in a particular (pt,v) bin. 



wi(2W3) 



,iOl+02-2<fe)\ _ 



(WlVl) (w 2 V 2 )vi(p T ,y)- 



(15) 



Note that there is no weight for the "differential" particle 
labelled 1, for which we do not perform any phase space 
average. With the previously derived values of (w 2 v 2 ) 
and (w±vi), one obtains the differential flow v\(pT,y)- 
This differential flow vi(pT,y) will be denoted by v[ in 
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this paper, while we keep the notation v\ for the inte- 
grated value only. 



III. IMPLEMENTATION 

In this section, we show how to analyze directed flow in 
practice using the three-particle correlation method. The 
method proposed here is a straightforward generalization 
of the one introduced in Ref. which involves the 

formalism of cumulants and generating functions. One 
may believe at first sight that this formalism is a useless 
complication here. But in fact, it provides an efficient 
and elegant solution to the following problems: 

• Taking into account azimuthal asymmetries in the 
detector acceptance, which always exist, even if the 
detector has full azimuthal coverage. 

• Eliminating correlations due to detector effects. 

• Dealing with the combinatorics, i.e., averaging over 
all possible triplets. 



We define the cumulants in Sec. [II A. The generating 



functions used in analyzin g inte grated and differential 



flow are introduced in Sec. Ill B. Then we give interpo- 



lation formulas which can be used to extract the r elevan t 
cumulants from these generating functions (Sec. Ill C| ). 



Finally, in Sec. HID, we derive the standard statistical 



If the detector is perfectly isotropic, the last term van- 
ishes and the cumulant reduces to the two-particle corre- 
lation. With a realistic detector, however, the cumulant 
isolates the physical correlation by subtracting the con- 
tribution of detector effects. Similarly, the cumulant of 
the three-particle correlation is defined as 



e 'i(0l+02-20 3 )\\ = / e i(0l+02-2tfe) 



J4>i \ / e i(02-20 3 ) 



J</> 2 \ / e iW>i-20 3 ) 



+2(e*)(e^ 2 )(e- 2 ^ 3 ). (17) 
As in the previous case, this isolates the genuine three- 
particle correlation from effects of detector inefficiencies, 
and from spurious correlations induced by detector ef- 
fects. 

If the acceptance is almost azimuthally symmetric, 
then all acceptance corrections are taken into account 
automatically by the cumulant expansion: in order to 
obtain the flow, one simply need replace the left-hand 
side (lhs) of Eq. (|) by the cumulant Eq. (0). If de- 
tector asymmetries are stronger, a multiplicative factor 
appears in the rhs of Eq. (||) relating the correlation to 
the flow. This correction is derived in appendix |A|. 



errors on both integrated and differential flow. 



A. Cumulants 



Even if the detector has full azimuthal coverage, its 
acceptance is not perfectly isotropic, so that averages 
like (e 1 ^ do not strictly vanish. A general way to take 
into account such effects consists in using cumulants. For 
instance, the cumulant associated with the two-particle 
correlation Eq. (||) is defined as 

(/ e i(^-M\ = (eH*t-**)\ - (e^ 1 ) (e-^ 2 ) . (16) 



B. Generating functions 

Generating functions provide an elegant way of sum- 
ming over all possible n-uplets in a given event. For a 
given event with M particles seen in the detector, we 
define the following real-valued function of two complex 
variables z\ = X\ + iy\ and z 2 — x 2 + iy 2 , 



M 

g(z U z 2 ) = n 

M 

= n 



i + !HM ( z * e ^ + zie -^- 



^ {z* 2 e 2 ^ + z 2 e- 2 ^ 



1 + ^jf- (2a:ioos(^) + 2j/i sin(^-)) + (2x 2 cos(20,) + 2y 2 sin(2^)) 



(18) 



where z\ = x\ — iyi and z\ = x 2 — iy 2 are the com- 
plex conjugates of z\ and z 2l respectively, and w\(j) and 
w 2 (j) are the weights mentioned in Sec. [IE. For sake 



of simplicity, we drop these weights from now on, unless 
otherwise stated. 
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The generating function G(z\, z 2 ) generalizes the gen- 
erating function G n (z) introduced in Ref. |fl7|| : the latter 
involved only one Fourier harmonic at a time, while we 
are now mixing two Fourier harmonics. More specifically, 
we recover the results of our earlier work in the limiting 
cases when either z\ or Z2 vanishes: G(z,0) = Gi(z), 
G(0,z) = G 2 (z). 

Neither the generating function G(z\, z 2 ), nor the com- 
plex numbers z%, z 2 have a physical meaning. They are 
a formal trick which allows us to extract azimuthal cor- 
relations to all orders, if necessary. This is done by av- 
eraging G(z\, z 2 ) over events (we denote this average by 
(G(zi, z 2 )}), and then expanding in power series of z%, 
z*, z 2 and zj. For instance, the coefficient of z* 2 ' z 2 in 
the expansion is 



coefficient of z* 2 z 2 in the expansion is 



(G( Zl ,z 2 )) = ■■■ + 



zfz 2 
M 3 



\j,k,i 



■20i) 



(19) 



The sum runs over nonequivalent triplets, i.e., j < k. 
The values of (j, k, I) are all different, so that auto- 
correlations are automatically avoided. Since there are 
M(M — 1)(M — 2)/2 nonequivalent triplets, one obtains 
for large M 



(G(z u z 2 )) = ■■■ + 



Z l Z 2 / e i(0 1 +0 2 -20 3 ) 



(20) 



One recognizes here the three-particle correlation Eq. (g) , 
averaged over triplets of particles and over events. 

In our averaging over events, we have assumed that 
the number of particles M is the same for all events. 
Although the method can accommodate small fluctua- 
tions of the multiplicity M, this is a possible source of 
error fTij] . Our recommendation is the following: in a 
given centrality bin, where the total number of detected 
particles M to t fluctuates from one event to the other, 
choose a fixed number of particles M < M tot to con- 
struct the generating function Eq. ([l8]). Although using 
only a fraction of the total number of particles results 
in a loss in statistics, this avoids the uncontrolled effects 
due to fluctuations in the multiplicity. 

Note that the generating function previously intro- 
duced in Refs. [jl7[ pfj|| was a function of only one complex 
variable z. Here, we need two independent variables Z\ 
and z 2 because we are mixing two different Fourier har- 
monics of the azimuthal distribution. 

Generating functions are not only a convenient way 
of summing over all n-uplets of particles. They also al- 
low one to construct easily cumulants of arbitrary order. 
The generating function of cumulants for integrated flow, 
C(zi,z 2 ), is defined by 0, || 

C(z 1 ,z 2 ) = M((G(z 1 ,z 2 )) 1/M -l) . (21) 

Expanding in power series of Z\ , z*, z 2 and z£ , one ob- 
tains cumulants of arbitrary order. In particular, the 



(C(z u z 2 )) = ■■■ + 



Zl z 2 



; »(01+02-20 3 )\\ + 



(22) 



An explicit calculation using Eqs. (O) and ([H]) shows 
that the expression 

f (e^i+^-^Jthus defined coin- 
cides with Eq. ( |l7j ) in the limit of large M. 

While the generating function C(zi, z 2 ) is real- valued, 
the cumulant ^(0i+<fe-2<fe)^ defined by Eq. @ is in 
general a complex number. However, the imaginary part 
results from detector effects and statistical fluctuations, 
and only the real part is relevant. For sake of brevity, 
we denote this cumulant of three-particle correlations by 
c{3} in the following: 



;{3} = Re 



+02-20 3 ) 



(23) 



where Re denotes the real part. Restoring the weights, 
this cumulant gives an estimate of the weighted inte- 
grated directed flow, which we denote by (iuivi{3}) [see 
Eq. (pi)]: 



c{3} = (wii>i{3}) 2 {w 2 v 2 ) 



(24) 



where the integrated elliptic flow (w 2 v 2 ) comes from an 
independent analysis. 

Let us now turn to differential flow. We shall denote 
by 1/) the azimuth of the differential particle under study, 
and by v[ its flow, v{ = (e^ - *")}, and call it a "proton" 
(although it can be any type of particle). In opposition, 
we call "pions" the particles used to estimate integrated 
flow. 

The overall procedure in the analysis is quite similar 
to the analysis of integrated flow. We first introduce a 
generating function of the azimuthal correlations between 
the proton and the pions. It is given by the average value 
over protons of e"^G(zi, z 2 ), where G(z\, z 2 ) is evaluated 
for the event where the protons belong: 



(e i +G{zx,z 2 )) = (e i +)+zx (e*^~^) 



(25) 



Note that in the averaging procedure, an event with two 
protons is counted twice, while an event with no proton 
does not contribute. 

Then we define a generating function of the cumulants 
for differential flow, V(zi, z 2 ), by [O, p0| 



T>{zi,z 2 ) 



(e^G( Zl ,z 2 )) 
(G{z u z 2 )) ■ 



(26) 



where, in the denominator, (G(zi,Z2)) denotes an aver- 
age over all events. The cumulants are the coefficients 
in the expansion in power series of this generating func- 
tion. As in the case of integrated flow, they are in general 
complex numbers, but only the real part is relevant from 
the physical point of view. For instance, the coefficient 
of z\z 2 defines the cumulant of order 3 which we shall 
use to extract differential flow: 



d{3} = Re (A 



3 i(V+0i-20 2 ) 



(27) 
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Note the similarity between this expression and Eq. (23). 
In fact, the cumulant d{3} shares the same features as 
c{3}: it is free from detector effects, and reflects physi- 
cal three-particle correlations, due either to direct three- 
body correlations, or to flow. Restoring the weights, 
the relation between the cumulant and flow reads [see 
Eq. @]: 



d{3} = (wivi) (W2V2) v[{3}. 



(28) 



We denote by ^{3} this estimate of v[ obtained from 
three-particle correlations. 



C. Interpolating the cumulants 

In this section, we show how to extract the cumu- 
lants for integrated and differential flows, Eqs. (^J) an d 
( p7| ) , numerically, from the computation of the generating 
function Eq. ( |ilf ) for various values of z\ and z 2 . We in- 
troduce the interpolation points (zi_ p ,Z2.q) — (xi,p+iyi,jj, 
X2,q + iy% q ) with 



xi, P = r cos y— J 

X2,q = fQCOS ^— J 



yi,p 

V2,q 



ro sm [—) 



tq sm 



for p — 0, . . . , 7 and q — 0, . . . , 3, and where ro is a 
real number, which must be neither too large, otherwise 
the error due to higher order terms in the power-series 
expansion of G(z\,Z2) is large, nor too small, to avoid 
numerical errors. 

To obtain the cumulants, one should for each event 
choose randomly M particles among the M to t detected, 
and, with the particle azimuths (f>j (and possibly with 
their transverse momenta and rapidities, if px and/or 
y-dependent weights are used) compute the generating 
function, Eq. (|l8|), at the points, Eq. (p9|). Then one 
must average the values G{z\ tPl Z2, q ) over events, and cal- 
culate the generating function of cumulants, Eq. (21). 



We denote by C VA the values of the generatin g fu nction 
C{z\,Z2) evaluated at the interpolation points (|29|): 



Cp,<J = C {Z\^ p , Z2,q) ■ 

From this quantities, we then build 

(Cp) x = (Cp,0 — Cp,2) 

{C p )y = (C Pi 3 — C P: i) 



(30) 



(31) 



which correspond to the real and imaginary parts of the 
derivative of C(zi,Z2) with respect to z 2 . Finally, the 
third order cumulant we are interested in, c{3}, is given 
by 

43} = [(Co), - (C\) v - (C 2 ) x + (C 3 ) y 

° + (C 4 ) x - (C 5 ) y - (C e ) x + (C 7 )y] . (32) 



Consider now differential flow. We denote by D Pi9 the 
values of the generating function 2?(z 1 ,z 2 ) evaluated at 
the interpolation points (C2T 



D ps = (D Pi q) x + i(D Ptq )y = V (zi iP , Z 2 ,q) , (33) 

where in fact one only need use even values of p: the in- 
terpolation of the "differential" cumulant, which is only a 
second order derivative, requires less points than the "in- 
tegrated" cumulant, which is a derivative of third order. 
With the quantities Eq. (|33|), we then build 

(D P ) X = j— [{D Pt0 ) x - (D Pi2 ) x + {D p ,i)y - (D Pj3 ) y ] , 
4 r Q 

{D p )y = - [(D p>0 )y - (D Pj2 ) y + (D Pt 3) x - (D Pt i) x ] , 
and the cumulant d{3} is finally given by 

43} = ^ l(D ) x - (Da), - (D 4 ) x + (D 6 ) y ] . (34) 

Naturally, one may prefer using a different interpola- 
tion scheme. In any case, one should check that the final 
values, c{3} and d{3}, do not depend on the parameters 
introduced in the interpolation: here, one should try dif- 
ferent values of r , and check the stability of the result. 



D. Statistical errors 

The standard deviation of the cumulant c{3}, Eq.(p3| 
is given by 



(Sc{3}f 



W) 2 (^ 2 ) 
2Af 3 7V cvts 



(2 + A X { + 2 X 2 2 + 4x1X2 + Xi) 



(35) 



In this equation, xi and X2 are the resolution param- 
eters appropriate for directed flow and elliptic flow, re- 
spectively: 



Xr, 



(w„V n ) 



M. 



(36) 



From E q. (p5| ), and the relation between c{3} and the 
flow, Eq. (|24D, one can calculate the statistical error on 
{w\Vi). If one neglects the statistical error on V2, one 
recovers the simple estimate (|ll|) with unit weights, in the 
limit where xi 1 and \i "C 1- A more careful estimate 
must take into account the error on (W2V2} in deriving 
the error on (wivi). Statistical errors on the cumulant 
are Gaussian. Since the relation between the cumulant 
and the integrated directed flow, Eq. (|24|), is quadratic 
rather than linear, the resulting error bars on (w\Vi) are 
asymmetric when the error is large, as discussed in detail 
in appendix D of Ref. jl7j . 

Similarly, the standard deviation of the differential cu- 
mulant d{3}, Eq. (p7j), is given by 



(<«{3}) 2 = 



2M 2 N' 



(I + X1XI + X2) 



(37) 
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where N' denotes the number of protons used in the dif- 
ferential flow analysis. Using the relation (|28| ) between 
this cumulant and the differential flow, ^{S}, we easily 
obtain the statistical error on i^{3}. Since the analysis 
is done in a narrow phase space bin, one may reasonably 
assume here that the error on c?{3} is dominated by the 
error on ^^{3}, which yields 



<5i^{3}/£ui{3} = \J MN cvts /N': errors scale like the in- 
verse square root of the number of particles involved, and 
the determination of the integrated flow v\ involves a to- 
tal number of MN cvts particles while the differential flow 
v[ involves N' particles. 



<K{3} 



i v/i + xfvi + xf 



Xi 



Xi 



(38) 



This result can be understood simply in two limiting 
cases. When both xi and Xi are large compared to 
unity, the reaction plane $i? can be reconstructed ac- 
curately. The error on (cos(V> — &r)) estimated with N' 
values of i/> is then l/y/2N', to which the error reduces for 
Xi > 1 , Xi > 1- in the opposite case xi < 1, Xi < 1, 
comparing Eqs. (B5h and (p7\), one finds for unit weights 



IV. HIGHER ORDERS 

In the previous section, we have studied the three- 
particle correlation, which is the lowest order nontrivial 
result obtained with the generalized generating function. 
Higher orders can also be derived in a straightforward 
way. Expanding the generating function of cumulants 
C(zi,Z2), defined in Eq. (|l]), in power series of z x , z*, 
Z2, and z%, yields cumulants of arbitrary order: 



C(z 1 ,z 2 ) = ^2 Z ^\ Z \ Z } Z f (^c' (01+,,,+01 ~ 0,+1 ^+*+2(^ + * + i+-+^ + k + i-^ J+fc+ j + i <j>i+h+i+m))\ (39) 

The only relevant cumulants are those which are invariant under a simultaneous shift of all azimuthal angles <f>j — > 
4>j + a. Other cumulants vanish except for statistical fluctuations and detector effects. According to Eq. (|l9|), this 
shift is equivalent to the change of variables z\ — > z\e~ la , z 2 — > Z2e~ 2la ■ The only terms in Eq. fl39| ) which are 
invariant under this transformation are those with j + 21 = k + 2m. Any of these cumulants can be used to extract the 
flow. Until now, we have explored only a few possibilities: the case considered in Sec. lTJ is (j,k,l,m) — (2,0,0, 1); 
the cumulants used in Ref. p7| to extract v\ and V2 are those with j = k, I = m — and those with j = k = 0,l = m 
(denoted by ci{2j} and C2{21}, respectively). 

We now derive the relations between cumulants of arbitrary orders and the flow coefficients v\ and v 2 - For this 
purpose, we compute the average value of G{z\, z 2 ) in the presence of flow. We first average for a given orientation of 
the reaction plane $>r. For an arbitrary particle, by definition of v n , we may write (e m *|$/?) = v n e m ^ R . Replacing 
in Eq. ([l9]), dropping the weights for simplicity, and neglecting all nonflow correlations, we obtain 



i®R _L -v ii c — _1_ -,*n, oli'&n. _i_ »i c — 2i$ f 



. z^ie'^'+ziwie +z 2 v 2 e—' i +z 2 v 2 e 
(G(Zi,Z2)\&r) = 1 + — 



M 

exp (z* lVl e l<s>R + z x v x e"***) exp {z* 2 v 2 e 2l * R + z 2 v 2 e- 2 ^ R ) . (40) 



The next step is to average over We make use of the 
following formula 

exp (z*e* + = e ~ iq \\z\) J * (2|z|) ' (41) 

where I q is the modified Bessel function of order q. Ap- 
plying this identity to each term of Eq. (^0|) and integrat- 
ing over <j>^, one obtains 
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(G( Zl ,z 2 )) = / (G(z 1 ,z 2 )\$r) 



d<f>R 
2ir 



+oo / *2 \<? 

E (^i^)/2 9 (2|z 1 K)/ 9 (2|z 2 h)(42) 



q— — oo 



zi\ 2 \zi\ 



In the limiting cases Z\ = or z 2 = 0, only the term 
q = contributes to the sum in the rhs, and we recover 
(G(zi,0)> = Jo(2|2iK), (G(0,«a)) =7o(2|^|« 2 ), already 
derived in Ref. 0. 



Expanding the generating function of cumulants, 
C(zi,z 2 ) ~ In (G(zi, z 2 )), in powers of zi, z*, z 2 , z 2 , 
and identifying with Eq. ([59|), one obtains the relations 
between the various cumulants and flow. As expected, 
the only nonvanishing terms are those which satisfy the 
condition j + 21 = k + 2m derived above, and the corre- 
sponding cumulants are proportional to (viy +k (v 2 ) l+m , 
with an integer multiplicative constant depending on the 
values of j,k,l,m. To order zl 2 z 2 z 2 , for instance, one 
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obtains 



e i(01+<fe-203+204-2<A 5 )\^ _ -(vi) 2 (v 2 ) 3 



(43) 



If V2 is measured independently, this equation yields an 
estimate of Vi from 5-particle correlations, which we de- 
note by t>i{5}. The statistical error on t>i{5} is not much 
larger than the error on the earlier three-particle esti- 
mate Wi{3} if V2 is large enough (more precisely, if the 
resolution parameter \2 ~ V2VM is larger than unity). 
This may be useful at RHIC, not at SPS where v-i is too 
small. Cumulants with j + k > 2 in Eq. (39) involve 
higher powers of v± and are useless in the situation we 
are interested in, namely small values of v\ . 



correlations between four (or more) particles. All in all, 
three-particle correlations seem to be the most appropri- 
ate way to measure v± when it is small, and especially if 
V2 is strong: near the balance energy, at SPS, RHIC, and 
the forthcoming experiments at LHC. 
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V. DISCUSSION 

In the previous sections, we have presented a new 
method for analyzing directed flow (vi ) , through the help 
of an independent measurement of elliptic flow (^2). It al- 
lows one to measure both integrated and differential flow, 
the value of integrated flow being used in the differential 
analysis. 

Our method relies on a study of £/iree-particle corre- 
lations. Unlike standard methods, based on two-particle 
correlations, it is not biased by two-particle nonflow cor- 
relations, which are an important bias at ultrarelativistic 
energies. This can be checked experimentally by study- 
ing vi near midrapidity: while standard two-particle es- 
timates «i{2} generally do not cross zero, especially if 
they are not corrected for momentum conservation [p.3| , 
our estimate t>i{3} should naturally vanish at midrapid- 
ity, as expected. 

In experiments where analyses of directed flow are al- 
ready available, it would be interesting to compare these 
standard estimates from two-particle correlations, «i{2}, 
with our new estimate from three-particle correlations 
«i{3}. If they are in agreement (within statistical er- 
ror bars), it is a good hint that they indeed coincide 
with the true directed flow. If they differ, it is instruc- 
tive to study the centrality dependence of the product 
M(«i{3} 2 — «i{2} 2 ), where M is the event multiplicity. 
If the difference between «i{2} and «i{3} is due to two- 
particle nonflow correlations, this product should be ap- 
proximately constant fl23| (remember that two-particle 
nonflow correlations scale as 1/M.) If the product differs 
significantly from a constant, then another explanation 
must be looked for: the difference between i>i{2} and 
«i{3} may be due to fluctuations, either fluctuations of 
the impact parameter within a given centrality class of 
events [^2|, or, more interestingly, physical fluctuations 
of the flow event-by-event. 

The price to pay for eliminating nonflow effects is an 
increase in statistical errors, compared to standard two- 
particle methods. However, this increase is moderate, a 
factor of 2 or less at SPS and RHIC. This new method is 
thus much less statistics-demanding than those based on 



APPENDIX A: ACCEPTANCE CORRECTIONS 

In the case of a realistic detector, with a nonisotropic 
acceptance, the relations between the cumulants and 
flow, Eqs. ( p4| ) and (p8|), no longer hold. In this ap- 
pendix, we derive the relevant modifications taking into 
account the detector acceptance. 

In what follows, we assume that the various classes 
of events analyzed, for instance, the different centrality 
bins, are determined with an independent detector (e.g., 
a ZDC), which has a full azimuthal coverage (at least 
approximately). This is meant to make sure that the 
centrality assigned to a given event is not strongly corre- 
lated to the orientation of its reaction plane, which would 
bias the sample of events used in the flow analysis. 

To describe a detector, we introduce its acceptance/ 
efficiency function A(j,<fi,pT,y)i which is the probabil- 
ity that a particle of type j with azimuth (/>, transverse 
momentum px , and rapidity y be detected [ [l7[ |2(| . Obvi- 
ously, A(j, 4>,pt, y) will vary from a detector to another, 
and a "perfect" detector corresponds to A(j, <p,pn y) = 1 
for every particle type in the whole phase space. In prac- 
tice, A(J,<f>,pT,y) is proportional to the number of hits 
in a (<P,Pt, y) bin: to obtain its shape for a given detec- 
tor and for each (pT,y) bin, one only has to count the 
number of hits in each (j) bin while reading the data to an- 
alyze them, and in the end divide by the maximal number 
encountered. The Fourier coefficients of the acceptance 
function are defined by 

A n (j,pT,y) = J 2 J A{] l( t>,p T ,y)e- m ^. (Ala) 

These differential coefficients can be integrated, with ap- 
propriate weighting and a sum over the various types of 
particles used for the flow analysis, so as to describe the 
"integrated" acceptance of the detector: 

^2 f w (j>PT,y) A n (j,p T ,y)dp T dy 

a n [w] = — j . (Alb) 

^2 w U'PT,y) A o(j,PT,y)dp T dy 
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Note our introducing the weights w(j,pj<,y), which are 
of course the same as in Eq. (|lq), either W\{j) or w 2 (j). 
In the following, we assume that both weights are equal, 
so that there is only one set of a n coefficients. If two dif- 
ferent weights are used — a weight w± which maximizes 
Vi and the weight w 2 which was used to derive the refer- 
ence v 2 — one should keep track of the two different sets 
of coefficients a n [wi], a n [w 2 ] in the calculation which we 
now sketch. 

To compute the contribution of flow to the cumulant 
c{3}, we follow the same procedure as in Rcf. [ p7[ . Wc 
first average the generating function Eq. (18) over events 
with the same azimuth of the reaction plane <E>/j; then 
we average over For simplicity, we neglect nonflow 
correlations in the derivation. Denoting by (x\$r) the 
average of a quantity x for fixed $r, the average value 
of e m * for a fixed $r is given by jl7j 



p#0 



p—n 



oXHe**«, (A2) 



where a* = a_ n is the complex conjugate of a n . Thus, a 
nonisotropic acceptance mixes the various flow harmon- 
ics. In the case of a perfect acceptance, Eqs. (Al) show 
that all coef ficients a n vanish, except a® = 1, and the 
identity (|a|) reads (e mcl> \$ R ) = v n e in ^ R , which follows 
in a straightforward way from Eq. ( pi). 

Inserting the average value, Eq. ( [A2[) for n = 1 and 
n = 2 in Eq. (|l8|), one obtains the generating function 
averaged over events with the same orientation of the re- 
action plane, (G(z\, z 2 )\<&r). The latter must then be 
averaged over then one computes the cumulants us- 
ing Eqs. ( pi] ) and J39|), keeping only the real part. In 
particular, the third cumulant reads 



c{3} = Re ((1 - |« 2 | 2 ) (1 - lail 2 )" 



(oi - a\a 2 ) 



\a 3 - aia 2 \ 



+ («4 - a\) (a* - a* 2 ) 



v\v 2 , (A3) 



instead of Eq. ( |24] ) , where we have assumed that only v\ 
and v 2 are nonvanishing. When the detector is perfect, 



one recovers Eq. (gj). Even if the detector is not perfect, 
but nevertheless does not have t oo b ad an acceptance, the 
factor in front of v\ v 2 in Eq. (A3) will remain close to 
1, since the correction terms are at least quadratic in the 
a n coefficients. Thus, Eq. (^) remains a good approxi- 
mation, except for detectors with a very bad azimuthal 
coverage. 

The contribution of flow to the "differential" cumulant 
d{3}, Eq. (p7|), can be calculated along the same lines. 
In that calculation, one may assume that integrated and 
differential flows are measured using two different detec- 
tors: e.g., a large acceptance detector for integrated flow, 
and a smaller one, but with better particle identification 
or Pt determination, for differential flow. We thus denote 
by A! '{j,ip,PT,y) the corresponding acceptance function 
and by A ' h (j, pt, y) its Fourier coefficients defined as in 
Eq. ( Ala ). The different ial ac ceptance coefficients a' k are 
then defined as in Eq. f JAlbj ), without the weights and 
the summation over j (since one usually measures the 
differential flow of identified particles) and with the inte- 
gration over pt and y restricted to the phase-space region 
under interest (typically, one integrates over pt or y, so 
as to obtain v[ as a function of y or px, respectively). 

The average value over protons in the numerator of Eq. 
( p6| ) is then computed in two steps, first averaging over 
protons detected in events with the same orientation of 
the reaction plane, then averaging over <&r. The denom- 
inator of Eq. (pf), (G(zi, z%) | is the same as in the 
calculation of c{3}. Finally, d{3} is given by the coeffi- 
cient of z\z 2 in the expansion in power series of T>(z\, z 2 ): 

d{3} = Re [(1 - |ai| 2 )(l - |a 2 | 2 ) + (oi - ata 2 )V 2 * 

+ |a 3 — a\a 2 \ 2 + (a* 2 — a* l 2 )(aA c — a\)a' 2 * v^viv-2 
+Re [(1 - |ai| 2 )(a 3 - aia 2 )a3* 

-I- (a* - a* 2 )(ai - a*a 2 )ai v' 2 v\. (A4) 

A nonisotropic acceptance will cause interference be- 
tween the various differential flow harmonics: the mea- 
surement of directed differential flow v[ is perturbed by 
the elliptic differential flow v' 2 . It is worth noting that as 
soon as the acceptance of the detector used for integrated 
flow is perfect, Eq. (|A4|) reduces to Eq. (ps|). 
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